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We analyze the response of a nanomechanical resonator to an external drive when it is also 
coupled to a single-electron transistor (SET). The interaction between the SET electrons and the 
mechanical resonator depends on the amplitude of the mechanical motion leading to a strongly 
non-linear response to the drive which is similar to that of a Duffing oscillator. We show that the 
average dynamics of the resonator is well-described by a simple effective model which incorporates 
damping and frequency renormalization terms which are amplitude dependent. We also find that 
for a certain range of parameters the system displays interesting bistable dynamics in which noise 
arising from charge fluctuations causes the resonator to switch slowly between different dynamical 
states. 
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I. INTRODUCTION 

Nanoelectromechanical systems (NEMS) in which the 
transport electrons in a mesoscopic conductor are cou- 
pled to a nanomechanical resonator 1 have been studied 
extensively over the last few years. Prominent examples 
include resonators coupled to tunnel junctions, 2 4 single- 
electron transistors (SETs) 5 " 8 or quantum dots. 9 " 11 
The electro-mechanical interaction in NEMS is typically 
rather weak and although there are important examples 
where non-linear coupling plays a significant role, 12-14 in 
many cases a linear description is sufficient. 

When the electro-mechanical coupling is weak, the 
effect of the resonator on the average current flowing 
through the conductor can provide an extremely sen- 
sitive measure of the mechanical displacement. How- 
ever, the electrons also act back on the nanomechan- 
ical resonator and in the weak-coupling limit their ef- 
fect on the resonator is typically analogous to a thermal 
bath. 1 ' 2 ' 15 Fluctuations in the current give rise to Gaus- 
sian fluctuations of the mechanical resonator character- 
ized by an effective temperature, which can be quite dif- 
ferent to the thermodynamic temperature of the system's 
surroundings. 6 " 9 

If the electro- mechanical coupling is increased then this 
simple picture inevitably breaks down and a more com- 
plex dynamics emerges. 15 " 18 The resonator dynamics can 
become highly non-linear, even if the electro-mechanical 
coupling itself remains linear, and the fluctuations can 
become no n- Gaussian though exactly what happens de- 
pends strongly on the type of charge transport involved. 
For example, for a resonator coupled to a normal state 
SET, the electrons always damp the mechanical motion 
on average, ' but the resonator probability distribu- 
tion nevertheless gradually changes from a Gaussian to 
a bimodal form as the coupling is increased. 16 ' 20 In con- 
trast, if the conductor tends to transfer energy to the 
mechanical resonator (as is the case for an appropri- 
ately tuned superconducting SET), increasing the cou- 
pling leads to a transition in the resonator dynamics 



which change from fluctuations about a fixed point to 
a state of self-sustaining oscillations 15 ' 17 ' 18 when the in- 
trinsic damping of the mechanical resonator is no longer 
sufficient to balance the energy transferred by the current 
flowing through the conductor. 

When the mechanical component of a NEMS with 
weak electro-mechanical coupling is driven to large ampli- 
tudes its influence on the charge dynamics of the conduc- 
tor is greatly enhanced. 21 ' 22 The resulting change in the 
charge transport also necessarily affects the way in which 
the charges act back on the mechanical system leading to 
a feedback process which can also generate strongly non- 
linear mechanical dynamics. Such effects have recently 
been investigated theoretically and seen experimentally 
in suspended carbon nanotube systems. 23 " 26 The effec- 
tive enhancement of electro-mechanical coupling in the 
presence of driving has been also been used in a novel 
form of force microscopy: when a cantilever which is 
capacitively coupled to quantum dots in a nearby sub- 
strate is strongly driven the resulting back-action on the 
mechanical dynamics can be used to infer information 
about the electronic structure of the dots. 

In this paper we use a simple model system consisting 
of a nanomechanical resonator linearly coupled to a nor- 
mal state SET 7 ' 16 ' 19 to explore how even very weak linear 
electro-mechanical coupling can give rise to a strongly 
non-linear response when the resonator is driven close to 
resonance. In the weak-coupling limit and in the absence 
of driving, the SET acts on the resonator like a thermal 
bath with an effective temperature proportional to the 
bias voltage; it also damps the mechanical motion and 
renormalizes the frequency of the resonator. 7 We find 
that for drives above a certain threshold the mechani- 
cal response as a function of frequency becomes strongly 
non-linear and the mechanical system displays many of 
the characteristics of the Duffing oscillator: frequency 
pulling, a strongly asymmetric line shape, hysteresis 
and bistability. Exploiting the fact that the underlying 
electro-mechanical coupling is very weak, we describe the 
effect of the SET on the resonator in terms of a simple 
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model which includes damping and frequency renormal- 
ization terms which are both amplitude dependent. 11 ' 28 
We find that a calculation of the average mechanical re- 
sponse as a function of drive frequency using these two 
quantities leads to results which are in very good agree- 
ment with a full Monte- Carlo simulation of the coupled 
dynamics. 

For a range of drive frequencies and amplitudes the av- 
erage mechanical response we calculate predicts the co- 
existence of high and low amplitude states. Monte-Carlo 
simulations show that in most cases the system spends 
all its time in just one of the two available states (which 
one depends on the initial state of the resonator). How- 
ever, we also find that there exists an interesting regime 
of bistability where the noise in the system is able to 
shift the system back and forth between the high and 
low amplitude states even when the two states are still 
quite different in terms of both the corresponding phases 
of the driven mechanical oscillations and average currents 
flowing through the SET. In this case the switching be- 
tween the two states is extremely slow compared to all 
the other time-scales of the system and could be detected 
in practice through a characteristic enhancement of the 
low frequency current noise. 

The structure of this paper is as follows. In Sec. II we 
outline our model for the driven SET-resonator system 
and derive master equations for the full coupled dynam- 
ics. Next, in Sec. Ill, we describe the behavior in the 
regime where the mechanical amplitude remains small 
and the resonator dynamics can be described using an ap- 
propriately modified version of the effective thermal bath 
description which applies to the undriven system. Then, 
in Sec. IV, we analyze the non-linear dynamics which oc- 
cur when the strength of the drive is increased. We calcu- 
late the effective damping and frequency renormalization 
of the resonator as a function of the amplitude and use 
these quantities to calculate the average mechanical re- 
sponse as a function of the drive frequency, comparing 
the results with those obtained from numerical simula- 
tions. In Sec. V we investigate a bistable regime where 
noise induced switching between two different states of 
the mechanical system occurs at a rate which is much 
slower than the internal dynamics of both the SET and 
the resonator. Finally, in Sec. VI, we present our conclu- 
sions. Brief details on the Monte- Carlo simulations are 
given in the Appendix. 



II. MODEL SYSTEM 

The SET-resonator system we consider is sketched 
schematically in Fig. 1. The SET island is coupled to 
the left and right leads by tunnel junctions with equal 
capacitances, Cj, and a bias voltage, V is assumed to be 
applied symmetrically. A gate electrode is used to tune 
the operating point of the island. The island capacitor 
consists of one plate which is mechanically compliant, 
giving rise to a position dependent capacitance, C g (x). 



v 
2 



v 

2 



gate 




FIG. 1. (Color online) Circuit diagram of the SET resonator 
system. The SET is coupled by tunnel junctions to two leads. 
A gate capacitor is used to tune the operating point, one 
plate of which is mechanically compliant to provide electro- 
mechanical coupling. 



This means that as the plate moves the operating point 
of the SET changes. Thus, as the charge on the island 
fluctuates, the electrostatic force on the plate changes, 
giving rise to electro-mechanical coupling. As long as 
the displacement of the resonator, x, is small compared 
to the distance, d, between the resonator and the SET 
island when the two are uncoupled, we can make a linear 
approximation for the dependence of the gate capacitance 
on the position of the resonator, 



C g {x) = C° (l - X - 



(1) 



where C® is the capacitance when the x = 0. 

The dynamics of the SET are determined by the rela- 
tive sizes of the energy scales in the system: the charging 
energy of the island, Ec = e 2 /2Cs (where is the to- 
tal capacitance of the island ), the thermal energy, k^T, 
and the energy scale of the bias voltage, eV. The Hamil- 
tonian for the system with n-charges on the SET island 
can be written as 



H n = Ec 



2nn) 



(>-§) 



P 

2m 



■xF(t), 
(2) 

where nP g = CgV g /e, p is the resonator momentum, ujo 
and m are the frequency and the mass of the resonator 
and F(t) is the external drive. We work in a regime where 
Ec ~ eV ^> /cbT, which means that only two charge 
states are accessible to the system, assuming < n° g < 1, 
these states are n = 0,1. 

There are four electron tunneling processes which can 
change the charge state of the island. We denote the 
rates for these processes by ^l(r) wnere +( — ) represents 
transitions which go in the same (opposite) direction to 
the applied bias, and L(R) denotes transitions at the left 
(right) junction as shown schematically in Fig. 1. The 
rates can be calculated within the orthodox model, 30 and 
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in the zero-temperature limit are given by 



L(R) 



AE L(R) 

e(^x«n)- 



J L(R)J Rje 2 



(3) 



where 0(-) is the Heaviside step function and Rj is the 
junction resistance. The AE terms are the free energy 
difference associated with each of the transitions, 



where, 



AE^ = ±E L ± moj^Tox, 
AE^ = ±E R =F mulxox, 



E L = Ec(l-2n g ) + ^- 
pV 

E R = -E c (l-2n° g ) + —, 



(4) 
(5) 



(6) 
(7) 



where {., .} is a Poisson bracket. Here we have defined the 
joint probability distributions ^o(i) v \ i) f° r the SET to 
have n = 0(1) and the resonator to have position, x, and 
velocity, u, at time t. The Hamiltonians #o(i) are the 
dimensionless form of Eq. (2) with n = 0, 1. 

The master equations [Eqs. (10) and (11)] together pro- 
vide a simple model for the SET-resonator system. The 
description is readily generalized to include the damping 
and thermal fluctuations of the mechanical resonator due 
to its interactions with its surroundings apart from the 
SET. However, we shall not include such effects explic- 
itly in our analysis and assume that the interaction with 
the SET electrons dominates the damping and fluctua- 
tions of the mechanical resonator. 8 ' 11 



III. LINEAR RESPONSE 



and we have introduced xq = —2n®Ec/mujQd, the change 
in the equilibrium displacement of the resonator when the 
SET changes between the charge states n = and n = 1. 
The bias voltage term, eV/2, accounts for the change in 
energy of the leads associated with the tunneling of an 
electron. 31 

We now adopt a dimensionless description: we intro- 
duce t = t/r, where r = (r+ + T+)- 1 = 2Rje 2 /eV, and 
the scaled position, x = x/xq. The dimensionless res- 
onator frequency is e = uoqt and the scaled tunnel rates 
are 

ff = B(±A L ± kx)(±A l ± kx), (8) 
f | = e(±A R =f kx)(±A r =F kx), , (9) 

where n = rmjj§x\j eV is the dimensionless electro- 
mechanical coupling strength and A L ^ = E L ^/eV. 
From now on we use the dimensionless forms of all quan- 
tities and so drop the tildes. 

The step functions in the expressions for the tun- 
nel rates have important consequences for large \x\. If 
x > x max , where x max = Ar/k, then the only possible 
allowed processes are those in which an electron tunnels 
onto the island r \ or so at most one tunneling event 
can occur until x < x max . The same kind of effect oc- 
curs when x < x m i n , with x m { n = —Al/k. In this case 
the only allowed processes are those in which an elec- 
tron tunnels off the island. Under such circumstances 
the resonator blocks transport through the SET, an ef- 
fect which is essentially the classical counterpart of the 
Frank-Condon blockade. 16 ' 20 ' 32 

The above expression for the Hamiltonians and tunnel 
rates can be used to write down classical master equa- 
tions which describe the evolution of probability distri- 
butions for the state of the SET and the resonator, 7 

P (x, v; t) = {H , P } + (T+ + r^A - (r+ + r+)p 0) 

(10) 

v, t) = {h u Pi} - (r+ + r^Px + (r+ + r+)p 0) 

(ii) 



When the dynamics of the system never take it into a 
region where the step functions in the tunnel rates, O(-), 
are important the system remains linear. This occurs as 
long as the condition k\x\ <C Al,r is satisfied for all of 
the phase space explored by the system. In practice these 
conditions are met provided both the electro-mechanical 
coupling is weak, k<1 and the drive is not too strong. 

Within the linear regime the master equations can be 
approximated as 

d l 2 ,mi^0 dP 
Po = [ex-f(t)}— -v — 

+ (A R - kx)P 1 - (A L + kx)P , (12) 

A-I^-')-/Wl^-^ 

- (A R - kx)Px + (A L + kx)P , (13) 

where f(t) = F(t)r 2 /mxo. Hence the full probability 
distribution of the resonator, P{x,v\t) = Po(x,v;t) + 
Pi(x,v;i), evolves according to 7 



P(x 1 v-t) = [e 2 x- f(t)] 



dP__ 2 dP L _ dP_ 
dv dv dx 



(14) 



These master equations can be used to find equations of 
motion for all moments of the position and velocity of 
the resonator, 



d_ 

dt 



IS 



(x n v m ) = // dxdvx n v m P(x,v-t). (15) 



We now analyze the behavior of the system for an ex- 
ternal driving force of the form, 



f(t) = f s'muj D t. 



(16) 



In the long-time limit the first moments of the mechanical 
resonator oscillate at the drive frequency, so we make the 
following ansatz for the fluctuating part of the position, 



5x — Ce 



(17) 
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where we have defined Sx = (x) — (x) ss with (x) ss the 
steady state of the undriven system with f(t) =0 (since 
the force has zero average when integrated over an integer 
number of periods). The equations of motion for the first 
moments of the system can be written as 



Sx = Sv, 
Sv = e 2 (SP-Sx) + f(t), 
SP = nSx - SP, 



(18) 
(19) 
(20) 



where Sv = (v) and SP = (Pi) - (Pi) ss with (Pi) = 
J J dxdvP\. Substituting the ansatz Eq. (17), into the 
expression for SP, and taking the Fourier transform we 
obtain 



SP(u) 



c S(u + uj d ) c J(uj-u d ) 



ILJ 



hence 



SP(t) = 



1 



-(Sx — Sv). 



(21) 



(22) 



D 



This allows us to write, 



Sx = 



1 



J D 



Sx 



e 2 K 



1 



■6v + f(t), (23) 



D 



which is simply the equation of a driven harmonic oscil- 
lator with renormalized frequency and damping due to 
the SET. These are given by, 



^eff 



1 



1 



• UJ 



D 



7eff 



• UJ 



(24) 



D 



These quantities take on a very similar form to those 
in the undriven system, 7 but now it is the frequency of 
the drive which enters these expressions in place of the 
natural frequency of the resonator. 

Solving Eq. (23) we obtain the coefficient from our 
ansatz, Eq. (17), 



C 



2(uj 2 ^-uj 2 d ) 2 ^2^uj 2 d • 



(25) 



The amplitude of the position oscillation of the resonator 
is then given by, 



2\C\. 



(26) 



To quantify the limits on the linear theory we introduce 
the critical amplitude, 



1 

2k' 



(27) 



As long as A x < A c the response as a function of fre- 
quency is a Lorentzian, centered around the renormalized 
frequency of the resonator. 

For an undriven resonator the fluctuations in the 
position and velocity are Gaussian in the weak cou- 
pling limit and can be described by invoking an effec- 
tive temperature. For the driven case a similar result 



is found, but the effective temperature is not a con- 
stant. The equations of motion for the variances (e.g 
Sx 2 = (x 2 ) — (x) 2 ) take the form, 



Sx 2 = 2Sxv, 
Sv 2 = 2e 2 (Svi — Sxv). 



(28) 
(29) 



This set of equations is completed by, 

Sxv = e 2 (Sx x - Sx 2 ) + Sv 2 , (30) 

Sxo — Svn — kSx 2 + ArSxi — AlSxo, (31) 

Sxi = Svi + kSx 2 — ArSxi + AlSxo, (32) 

= — e 2 Sxo — kSxv + ArSvi — AlSvq — e 2 (Pi)(Po), 

(33) 

5i>! = -e 2 8xi + k5xv - A R dv! + A L Sv + e 2 (P 1 )(P }, 

(34) 

where (Pq) = 1 — (Pi) and we have defined, for example, 
Sxo = J J dxdvPox — (x)(Po). Although these equations 
can be solved analytically we do not give the solution 
here as it is rather cumbersome. 

The dynamics of the variances is simplest in the adia- 
batic limit where ujd <C 7 e ff- In this case the resonator 
can be thought of as relaxing to a thermal distribution 
at each point during the drive cycle. In this case the 
fluctuations in the position and velocity of the resonator 
can be described in the same way as the undriven case 
by an effective thermal distribution with 



eV 



(PiWXPoW), 



(35) 



though now the term (Pi(t))(Po(t)) oscillates according 
to Eqs. (21) and (25). The oscillations in the SET charge 
therefore generate oscillating components in the effective 
temperature (and hence the variances Sx 2 and Sv 2 ) at ujd 
and 2cj£>. When ujd ~ 7eff the oscillations get smaller 
as the resonator cannot relax fast enough to follow the 
oscillations in the effective temperature. If the drive is 
much faster than then the oscillations are washed out 
and the variances of the resonator are set by the effective 
temperature averaged over the drive period. 

The analytical results in the linear regime can be com- 
pared with Monte- Carlo simulations (see the appendix 
for more details) of the system dynamics. 35 An exam- 
ple of the full position distribution over one period of the 
drive is shown in Fig. 2(a). As expected, the mean is well 
described by the linear analytic result from Eq. (17). In 
Fig. 2(b) we compare the numerically simulated variance 
in the distribution over one drive period to the solution of 
Eqs. (28) — (34). The two curves show good agreement, 
the small deviations occur because a few of the numer- 
ical trajectories enter the non-linear region. In Fig. 2, 
we have chosen uoD/left ~ 1-4, and so the oscillations 
in the variance are visible, but are smaller than the am- 
plitude of the oscillations in the corresponding effective 
temperature in the adiabatic limit [Eq. (35)]. 
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FIG. 2. (Color online) Oscillations in both the mean and 
variance of the position distribution over a period of the drive. 
The parameters are e = 0.3, n = 0.05, fo = 0.05 and ujd = 
0.0027T. In (a) we show the full position distribution (obtained 
from a Monte- Carlo calculation) and the analytic solution for 
the mean as the dashed (white) line. In (b) we show the 
variance, the solid line is the numerical result, the dashed 
line shows the solution to Eqs. (28)- (34). Note that for all 
numerical calculations we chose Al, Ar so that the system 
is at the charge degeneracy point, (Pi) = 1/2, where the 
oscillations in (Pi)(Po) at frequency ujd vanish. 




IV. NON-LINEAR RESPONSE 

Even if the electro-mechanical coupling is extremely 
weak, k « 1, a slight increase in the drive strength can 
be enough to reach a regime in which the amplitude of 
the oscillations exceeds A c . Thus, the resonator explores 
regions of phase space where the effects of the Frank- 
Condon blockade are important. When this happens the 
tunneling processes in the SET start to become signifi- 
cantly modified by the step function constraints, the lin- 
ear theory breaks down and a richer, more complex dy- 
namics emerges. In this section we will examine the re- 
sponse of the resonator to a drive which is strong enough 
to lead to non-linear behavior and develop techniques to 
describe the dynamics in this regime. We start with a 
very simple approximate way of including the effects of 
the step functions. This approach provides an intuitive 
description of the non-linear dynamics of the resonator 
which is qualitatively correct. A more detailed calcu- 
lation of the amplitude dependence of the damping and 
frequency is then presented which is in good quantitative 
agreement with Monte- Carlo simulations. 



A. Reduced coupling approach 

When the system is outside of the linear region, the 
effective damping and frequency shift of the resonator 
will be modified and the expressions in Eq. (24) are no 
longer valid. The main effect of the step functions in the 
tunnel rates is to block electrons from traveling through 
the SET and when electrons cease to flow their back- 
action on the resonator will stop too. Assuming that 
all electron tunneling stops whenever x > x max or x < 
^min, ' we can account for the consequent reduction in 
the damping and frequency shift simply by weighing the 
effective coupling strength by the proportion of time the 



FIG. 3. (Color online) Amplitude dependent (a) damping and 
(b) frequency of the resonator. Solid (black) lines show the 
solution using Eq. (50) dashed (red) lines show the results 
using Eq. (39). The parameters used are e = 0.3, k = 0.05, 
fo = 0.015 and oj d = 0.29. 

resonator spends inside this region. 35 

We start by assuming that the position of the resonator 
varies harmonically at the frequency of the drive 36 , 

Sx(t) = A x sin(cjc>t). (36) 

Using this position dependence and the assumption that 
the influence of the SET on the resonator is switched off 
when x > x max or x < x m i n , we can define an amplitude 
dependent effective coupling, 

27r Jo 

(37) 

which is readily integrated to give 

f — arcsin ( tt 1 *-) A x > ±- . . 

The renormalized coupling gives rise to an implicit am- 
plitude dependence in the damping and frequency shift, 

(39) 

Examples of these functions are shown as dashed lines 
in Fig. 3. In the region where the amplitude is still in the 
linear region, A x < A C1 the damping and frequency shift 
retain their linear values. Above A c we see the damping 
decrease towards zero and the frequency move towards 
the bare frequency, e. In the large amplitude limit the 
system spends less and less time inside the linear region 



6 



where the SET electrons damp the motion. We note 
that at sufficiently large amplitudes the damping of the 
resonator due to sources other than the SET electrons 
(which we neglected) will necessarily stabilize the dynam- 
ics even if, as we assume here, such damping is very small 
compared to 7 e ff (0). 

Solutions for the amplitude of the oscillations are ob- 
tained self-consistently from the expression 



At 



/o 



(40) 



which is derived from the equation of motion for a driven 
resonator with amplitude dependent damping and fre- 
quency. Equation (40) can have either one or three so- 
lutions depending on the choice of parameters. To find 
which of these solutions are stable we derive an equa- 
tion of motion for A x . 28 This is given by, 



A x 



1 

A Ti 



SxSx 



SvSv 



D 



(41) 



The problem is greatly simplified if we assume A x varies 
much more slowly than the drive force (to be expected 
when the coupling is weak). We therefore introduce a 
period averaged amplitude, 28 



A.f _ / A-rn (ljty , 

which obeys the equation of motion, 



dA x 



g{A x ), 



(42) 



(43) 



where, 



9(A) 



7eff(A) 



/o 2 



A[(w* eS (A)-w%) + >&(A)w%] 



(44) 

The fixed point solutions are given by g(Af p ) = which 
is just Eq. (40). 

The fixed point amplitudes are stable when 28 ' 37 



dg 



dA r 



< 0. 



(45) 



Hence we find that where Eq. (40) has three solutions: 
the small amplitude solution is in the linear region (and 
hence stable), the intermediate amplitude solution is un- 
stable and the large amplitude solution is a new stable 
fixed point for the system (where the damping and fre- 
quency shift are reduced from the linear values). The 
presence of more than one stable amplitude is com- 
mon in driven non-linear systems, such as the Duffing 
oscillator. 37 However, in contrast to the standard Duffing 
oscillator 37 where the nonlinearity arises from the poten- 
tial (and the damping is always linear), the interaction 
with the SET charge leads to both non-linear damping 
and frequency terms. 



Using the non-linear damping and frequency shift, 
along with the conditions on stability of the solutions to 
Eq. (40), we obtain the response of the system, A x , as a 
function of drive frequency. An example of the resulting 
curve can be seen in Fig. 4. At low frequencies (below the 
shaded region labeled (a) in Fig. 4) and at high frequen- 
cies ujd > e, the system remains in the linear regime; the 
response is below A c . However, the response to the drive 
becomes stronger closer to resonance. In the shaded re- 
gion labeled (a) in Fig. 4, the amplitude grows beyond A c , 
and so the linear and non-linear calculations give differ- 
ent results. The linear calculation leads to a Lorentzian 
peak centered around cj e ff(0). However, in the non-linear 
case, the frequency shift becomes smaller (leading to a 
larger effective frequency) for A x > A c . This means that 
the drive frequency is further from resonance than in the 
linear case, and hence the amplitude is smaller in shaded 
region (a). 

In the shaded region of Fig. 4 labeled (b), the drive fre- 
quencies are close to (but below) e and a high amplitude 
solution exists because of a positive feedback mechanism. 
In this regime, as the amplitude grows beyond A c , the 
enhancement in the effective frequency brings the system 
closer to resonance with the drive (and the damping de- 
creases as the amplitude grows). However, the system 
does eventually stabilize when the amplitude becomes 
large enough that the system starts to move away from 
resonance, when cj e s(A x ) starts to increase beyond ujd- 

For drive frequencies larger than e, the high amplitude 
solution no longer exists. To see why, consider starting in 
the high amplitude state and then increasing ud so it is 
slightly larger than e. If the damping and frequency shift 
were not amplitude dependent the system would simply 
oscillate at a slightly lower A x , due to being more de- 
tuned from resonance. However, as the amplitude drops, 
the effective frequency shifts further from resonance (cj e ff 
decreases with decreasing amplitude) and the effective 
damping increases, reducing the amplitude of the oscilla- 
tions further: the system spirals back down to the linear 
branch. 



B. Amplitude dependent damping and frequency 

The simple argument presented above is able to give 
qualitative agreement with the numerical simulations, 
but to gain quantitative agreement we need to use a more 
accurate method to calculate the effective damping rate 
and renormalized frequency. 11 ' 28 We assume that a con- 
stant amplitude solution to the full non-linear expression 
exists, and use this to map the equations onto those of 
damped, driven harmonic oscillator. We can then iden- 
tify the terms which correspond to an amplitude depen- 
dent damping and frequency shift. 

We begin with the set of coupled equations for the 
first moments, including the non-linear position depen- 
dence of the tunnel rates (which enters through the step 
functions), but decoupling the second moments which 
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FIG. 4. (Color online) The solid (black) lines show the stable 
solutions to the self-consistency expression, Eq. (40), as a 
function of drive frequency, for fo = 0.02, e = 0.3 and k, = 
0.05. We also show the linear result as the dashed (black) 
curve and the amplitude at which the linear theory fails, A c , 
as the horizontal dot-dashed (blue) line. The shaded regions 
(a) and (b) are discussed in the text. 

are present, so that, for example, we assume (xi) = 
(x)(Pi). 28 We also approximate averages of a function 
of position, (/(#)), by the equivalent function of the av- 
erage f((x)). These approximations have been found to 
work well for a range of similar systems. 11 ' 28 ' 38 We there- 
fore obtain 

(i)=e 2 ((P 1 )-{x)) + f(t), (46) 
(A) = r+((s»(l - (Fx)) + r+((x))(P 1 ). (47) 

Back-tunneling (i.e. processes described by F~£( R \) are 
neglected as their contribution remains very small even 
within the non-linear regime. 35 

We now solve for (Pi(t)), by integrating Eq. (47) nu- 
merically using the same ansatz for (x(t)) as Eq. (36). 
Making a comparison between Eq. (46) and the driven 
oscillator equation we identify 

- e 2 (P 1 ) = 7eff (x) + A^ 2 ff (x) , (48) 

where we have defined the frequency shift, 

Aw e 2 ff = cv 2 eS - e 2 . (49) 

Multiplying Eq. (48) by either cos ojjjt or sincj^t, and 
integrating over one period of the drive we obtain, 18 

£ 2 p2<k/uj d 

7eff = t / Pi(t) cos UDtdt, (50a) 

( I>2tt/u} D \ 
l--|y Pi(t)smuj D tdt\ . (50b) 

These expressions correspond to the usual physical in- 
terpretation of the damping and frequency shift arising 
from the SET electrons: the damping is due to the out 
of phase (with respect to x) component of the average 
island charge, while the frequency shift is due to the in 
phase component. 11 ' 18 ' 20 



In Fig. 3 we plot the amplitude dependence of the 
damping and frequency shift calculated using Eq. (50), 
along with the simple results calculated previously using 
Eq. (39). We see that Eq. (50) is always larger than the 
simple estimates. This is because the damping does not 
simply switch off as soon as the mean resonator position 
travels through the non-linear boundary as some electron 
tunneling processes can still occur even when (x) < x m i n 
or (x) > x max . 

The amplitude dependent damping and frequency shift 
can be used in the self-consistency expression, Eq. (40), 
to find the amplitude of the stable solutions. We now 
compare the predicted stable amplitudes to Monte-Carlo 
results for parameters which enter the non-linear region. 
Since the system has two stable solutions with different 
amplitudes the dynamics are more complicated than in 
the linear case. Figure 5 compares the amplitudes ob- 
tained using the effective damping and frequency shift 
with results obtained from Monte-Carlo simulations. In 
general there is very good agreement, though there is a 
small region where the low amplitude solution predicted 
from the amplitude dependent damping and frequency is 
not seen in the numerics. This small discrepancy is a re- 
sult of the fluctuations in the system which are included 
in the Monte-Carlo calculation. 

The numerical results are obtained by either increas- 
ing or decreasing the frequency progressively and in each 
case taking the state of the system after a long trajec- 
tory at a given frequency as the initial condition for the 
next frequency value. Calculated in this way the results 
show clear hysteresis: when the frequency is swept for- 
wards (from lower to higher values) the resonator always 
follows the high amplitude branch. However, when fre- 
quency is swept downwards (from higher to lower fre- 
quency) the system remains in the low amplitude state 
for a long time because although a high amplitude so- 
lution exists it is too far away from the low amplitude 
one to be reached by fluctuations on the time-scale of 
the simulations. However, as the drive frequency is pro- 
gressively reduced the amplitude of the high amplitude 
oscillations decreases and that of the low amplitude oscil- 
lations increases. Eventually the fluctuations in the low 
amplitude state are enough to take the resonator out of 
the linear regime and the low amplitude state becomes 
completely unstable and is no longer seen. 

Figure 5(b) shows the current (averaged over a large in- 
teger number of periods) as a function of drive frequency 
for both linear and non-linear cases obtained from Monte- 
Carlo simulations. For the linear case the current is sup- 
pressed in the region where the resonator is resonantly 
driven; as the amplitude of the oscillations increases one 
of the tunnel rates is suppressed leading to an overall 
reduction in current (the system is tuned to a gate volt- 
age which for the undriven steady state position corre- 
sponds to a current maximum, the charge degeneracy 
point). The current flowing through the SET reflects the 
non-linear behavior of the resonator and there is a clear 
dependence on the direction of sweep. For uod — ^ the 
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FIG. 6. (Color online) Variance of the x distribution averaged 
over an integer number of drive periods, (5x 2 ), for various 
drive strengths, the solid (black) curve is fo = 0.01 (linear), 
the dashed (purple) curve is fo = 0.015 (bistable) and the 
dot-dash (green) curve is fo = 0.02 (non- linear). These re- 
sults were obtained from Monte- Carlo simulations in which 
the frequency was swept downwards 



FIG. 5. (Color online) Amplitude of the resonator (a) and 
average current (b) as a function of drive frequency. Results 
from a Monte- Carlo simulation with a forward (solid (red) 
line) and backward (dot-dash (green) curve) frequency sweep 
for fo = 0.025 are shown together with the fully linear re- 
sults for fo — 0.01 (solid (black) curve). Notice that for each 
sweep the initial conditions for a given frequency were cho- 
sen to match the final state found for the previous value of 
the frequency. The results obtained using amplitude depen- 
dent effective damping and frequency shift [Eq. (50)] are also 
shown in (a) plotted as a (dashed (blue) curve). The other 
parameters used are e = 0.3 and k = 0.05. 

current is reduced almost to zero in the forward frequency 
sweep: the amplitude of oscillation of the resonator is so 
large that only one tunnel event occurs each time the 
resonators passes through the region where transport is 
allowed. 



V. BISTABILITY 

In contrast to the calculations of the effective damping 
and frequency shift the Monte- Carlo simulations capture 
the fluctuations in the resonator's dynamics. For the pa- 
rameters used in Fig. 5 fluctuations take the system from 
a low to a high amplitude state at a particular drive fre- 
quency, but the system is never able to switch back and 
forth between states of different amplitude. However, 
it is possible to find parameters where bistability does 
occur by adjusting the drive amplitude and frequency 
so that the 'low' and 'high' amplitude states are close 
enough that fluctuations can carry the system between 
the states in both directions. Bistability is common to 
a wide range of damped, driven non-linear systems, such 
as in the Duffing oscillator, but is also seen in other 
non-linear NEMS such as a resonator coupled to a su- 
perconducting SET. 18 ' 28 



The fluctuations in the state of the resonator provide 
important information about its dynamics. The variance 
of the position averaged over a large integer number of 
drive periods, (5x 2 ), is shown in Fig. 6 for different drive 
strengths. At low drives the dynamics are linear and the 
variance is approximately constant, with a shallow dip 
around the resonance. 39 The fluctuations are typically 
much larger for higher drive strengths when the system 
enters the non-linear regime as can be seen from the dot- 
dashed curve in Fig. 6 which corresponds to the sweep 
from high to low frequency in Fig. 5. However, by far 
the largest fluctuations are seen for the dashed curve 
in Fig. 6 which corresponds to an intermediate drive 
strength. A sharp peak in the variance suggests bista- 
bility, since a bimodal probability distribution typically 
gives rise to a very large variance. Figure 7 shows the 
average amplitude for the same parameters which give 
rise to the sharp peak in Fig. 6. For this drive strength, 
in contrast to the behavior seen in Fig. 5(a), there is a 
smooth crossover, independent of the direction of the fre- 
quency sweep, between the two different states predicted 
by the effective damping and frequency calculation. This 
smooth crossover in amplitude at co>£>/e ~ 0.98 and the 
corresponding peak in the fluctuations strongly suggest 
that the system has a bistable region and we can con- 
firm that this is indeed the case by examining individual 
Monte-Carlo trajectories. 

Examples of the dynamics during a single trajectory 
are shown in Fig. 8, where we show a trace of both the 
current through the SET and the resonator's phase. 40 To 
make the behavior clear each point in the trajectories is 
averaged over a timescale which is long compared to the 
drive period, but short compared to the rate at which 
the system moves between the two different states. It is 
clear that the dynamics are well described by a two state 
model, the system switches between two distinct values 
for both the current and the phase as a function of time. 
This allows us to split the trajectory into regions which 
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FIG. 7. Response of the system as a function of driving 
frequency, all parameters are the same as in Fig. 4, except 
fo — 0.015. The solid line shows the numerical response, the 
dashed line shows the analytically predicted stable solutions 
using the damping and frequency shift from Eq. (50). The 
system shows a bistability in the region where the linear so- 
lution reappears. 




FIG. 8. Example trajectories of the current through the SET 
and phase of the resonator as a function of time. The dashed 
lines show the point at which we choose to split the data 
into the high and low amplitude states. The parameters are 
e = 0.3, k = 0.05 and fo = 0.016. The drive frequency 
is chosen such that the occupation probabilities of the two 
states are equal. Each point in the plots is averaged over 200 
periods of the drive. 



correspond to the 'high' and 'low' amplitude oscillations 
of the resonator. The phase and current trace are well 
correlated, when the system is in the 'high' amplitude 
oscillation state the current is small and the phase is 
large, the opposite is true for the 'low' amplitude state. 

The transition rates between the two bistable states, 
Thi and T^, where h(l) labels the high (low) amplitude 
state, are given by 



hi 



1 

Th 



ih 



1 

n 



(51) 



where is the average amount of time spent in the 
state h(l). These rates 41 are readily extracted from the 
Monte- Carlo trajectories. Figure 9(a) shows the behav- 
ior of T h i(i h ^ as the drive frequency is swept through the 
bistability. The rates are much slower than all other 
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FIG. 9. (Color online) (a) Transition rates between the two 
states as a function of drive frequency, close to the bistability. 
The (blue) dots denote Yui the rate from the high ampli- 
tude state to the low amplitude state, the (red) crosses are 
T^. The insets show the resonator probability distributions 
at the drive frequencies labeled (i), (ii) and (iii) in the main 
panel. The parameters at point (ii) are the same as in Fig. 
8. (b) Corresponding zero- frequency Fano factor, calculated 
from the two state model, as described in the text. The points 
labeled (i), (ii) and (iii) match those in (a). 



timescales in the problem, for example 7 e ff ~ 0.005 for 
the low amplitude state and 7 e ff ~ 0.002 for the high 
amplitude state. This means that the system has time 
to relax in each of the metastable states and a two state 
model should be a good approximation to the dynam- 
ics on intermediate time-scales. 17 ' 42 The insets to Fig. 9 
show probability distributions for the resonator at the 
three marked drive frequencies. These are obtained stro- 
boscopically by taking points at a particular drive phase. 
We see that when Thi < Tih, illustrated in inset (i), 
the majority of the distribution is in the high amplitude 
state, which has a wide distribution in phase. At frequen- 
cies above the bistability [inset (iii)] the distribution is 
dominated by the 'low' amplitude state (which is essen- 
tially the linear solution here), since T^i > Fih- Between 
these two limits [inset (ii)] the two rates are approxi- 
mately equal, Fhi ~ T^, and the distribution contains 
significant contributions from both states. 

The noise in the current flowing through NEMS is 
known to provide important information about the dy- 
namics of the system. 17 ' 38 ' 42 ' 43 Assuming that the current 
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is measured on a timescale slower than the drive then the 
relevant quantity is the period averaged current 44 



Kt')dt'. 



(52) 



It is the current defined in this way which shows the 
bistable behavior seen in Fig. 8. Close to zero frequency 
the noise in this current will be dominated by the switch- 
ing rate between the two states, since these are the slow- 
est timescales in the dynamics. 17 ' 38 ' 42 The noise can the 
be quantified by the zero-frequency Fano factor of the 
current, 



F(0) 



(53) 



Assuming the system is well-described by a two-state 
model it is possible to find a simple analytic expression 

for F(0) 17 ' 38 ' 42 



F(0) 



2P h p l ({i h )-{i l )y 

T h i{I h )+Ti h (Ii) 



(54) 



where Pr{1) an d (Ih(i)) are the occupation probability and 
average current of the high (low) amplitude state respec- 
tively. These quantities can be calculated in a straight- 
forward manner from traces like those in Fig. 8. The 
predictions for the behavior of the Fano factor based on 
the two-state description are shown in Fig. 9(b). As is 
typical for this kind of dynamics, 17 ' 38 ' 42 ' 43 we find a large 
enhancement in the noise close to the bistability: switch- 
ing between the two metastable states causes large fluc- 
tuations in the current. The maximum in the Fano factor 
occurs for a drive frequency slightly larger than that for 
which the occupation probabilities of the two states are 
equal, (the point (ii) in Fig. 9 where Thi = ^ih)- As the 
frequency is increased the difference between the current 
in the two states, \(Ih) — (Ii)\ increases, as can be seen in 
Fig. 5(b), this contribution shifts the maximum in F(0) 
slightly away from the point where Ph = Pi. 



VI. CONCLUSIONS 

We have used a simple model system consisting of 
a mechanical resonator linearly coupled to a normal 
state SET to explore the non-linear dynamics which can 
arise in driven NEMS. Despite weak (linear) electro- 
mechanical coupling the resonator's dynamics become 
non-linear when it is driven sufficiently strongly. At large 
enough amplitudes the effect of the resonator on the SET 
charge dynamics can no longer be accounted for by a lin- 
ear correction to the tunnel rates and the charge trans- 
port is strongly modified. The modified charge dynamics 
leads to changes in the damping and frequency shift in- 
duced by the SET on the resonator leading in general 
to an amplitude dependence of these quantities. Such 



an amplitude dependence is generic in non-linear oscilla- 
tors and leads to the familiar phenomena of asymmetric 
frequency response, hysteresis and bistability. 

We have focused on a simple model in which the effects 
of finite temperature on the electron tunneling and the 
effects of the resonator's surroundings beyond the SET 
electrons are ignored. Whilst the model can be adapted 
to include all of these and other complicating factors it is 
nevertheless useful to work with a simplified description 
as it lays bare the mechanisms which give rise to non- 
linear behavior. 

For weak, linear, electro-mechanical coupling in NEMS 
one generally expects the electrical transport to act on 
the resonator like an effective thermal bath, giving rise 
to damping of the mechanical motion and Gaussian 
fluctuations. 2 ' 7 ' 8 ' 15 The cases where coupling between the 
electrons and the resonator instead gives rise to negative 
damping have become established as well-known excep- 
tions to this paradigm. 15 ' 17 ' 18 The way in which the ef- 
fective thermal description breaks down as the electro- 
mechanical coupling is increased has also been investi- 
gated carefully. 16 ' 20 Our work has explored a third case: 
in general one expects driving of the resonator in a NEMS 
device to lead (via the interaction with the transport elec- 
trons) to non-linear dynamics which also fall well outside 
the effective thermal bath description. 26 
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APPENDIX: MONTE-CARLO SIMULATIONS 

In this appendix we briefly summarize the numerical 
techniques used to obtain the Monte-Carlo results dis- 
cussed in the main body of the article. The basic idea 
is to simulate the dynamics obtained from Eq. (2) via 
Hamilton's equations, 



-e\x-n) + f(t). 



(55) 
(56) 



When electron tunneling is taken into account these 
equations are essentially stochastic, because the island 
charge, n, fluctuates between zero and unity according 
to the tunnel rates given by Eqs. (8) and (9). 

To perform Monte Carlo simulations of the trajectory 
of the system governed by these equations we evolve time 
in discrete steps of length At. At each time step, a 
uniform random number, r G [0,1], is generated and 
compared with the appropriate rates to see whether the 
charge state should be updated. For example, if n = 
and r is in the interval [0, (rj + T^)At] then the state 
of the SET island is changed n = 1, otherwise no change 
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is made. The resonator is then evolved over the timestep 
using Eqs. 55 and 56. 

The results shown in the main body of the article are 
obtained by appropriate averaging over single long trajec- 
tories. Starting from an arbitrary set of initial conditions 



the trajectory is evolved a time, to, long enough that the 
system has lost all memory of the initial condition (this 
is tested by ensuring that the results are not sensitive 
to changes in to). We then evolve the system along the 
trajectory recording running averages as required. 
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